/*
 * Reaction.h
 *
 *  Created on: 23 May 2011
 *      Author: allan
 */

#ifndef REACTION_H_
#define REACTION_H_

// C++ includes
#include <iostream>
#include <memory>
using namespace std;

// Eigen includes
#include <Eigen/Core>
using namespace Eigen;

// GeoReact includes
#include <Core/EquilibriumConstant.h>
#include <Core/Kinetics.h>
#include <Core/Multiphase.h>
#include <Utils/Types.h>

class ReactionEquation
{
public:
	ReactionEquation();

	ReactionEquation(const vector<string>& species, const vector<double>& stoichiometries);
	
	friend class Reaction;

private:
	/// The species that participate in the reaction
	vector<string> species;

	/// The stoichiometries of the species
	vector<double> stoichiometries;
};

class Reaction
{
public:
	Reaction();
	
	Reaction(const ReactionEquation& equation,
		const Multiphase& multiphase);
	
	Reaction(const ReactionEquation& equation,
		const EquilibriumConstant& K,
		const Multiphase& multiphase);
	
	template<class KineticsType>
	Reaction(const ReactionEquation& equation,
		const EquilibriumConstant& K,
		const KineticsType& kinetics,
		const Multiphase& multiphase);
	
	void
	SetEquilibriumConstant(const EquilibriumConstant& K);
	
	template<class KineticsType> void
	SetKinetics(const KineticsType& kinetics);
	
	const vector<string>&
	GetSpecies() const;
	
	const vector<double>&
	GetStoichiometries() const;
	
	const vector<Index>&
	GetSpeciesIndexes() const;
	
	const EquilibriumConstant&
	GetEquilibriumConstant() const;
	
	const Kinetics&
	GetKinetics() const;
	
	const double
	Stoichiometry(const string& species) const;
	
	const double
	EquilibriumConst(double T, double P) const;
	
	const double
	ReactionQuotient(const VectorXd& a) const;
	
	const double
	KineticRate(double T, double P, const VectorXd& n, const VectorXd& a) const;

private:
	/// The species that participate in the reaction
	vector<string> species;

	/// The stoichiometries of the species
	vector<double> stoichiometries;

	/// The indexes of the species in the reaction
	vector<Index> indexes_species;

	/// The equilibrium constant of the reaction
	EquilibriumConstant K;

	/// The kinetics of the reaction
	std::shared_ptr<Kinetics> kinetics;
};

template<class KineticsType>
Reaction::Reaction(const ReactionEquation& equation, const EquilibriumConstant& K, const KineticsType& kinetics, const Multiphase& multiphase) :
species(equation.species), stoichiometries(equation.stoichiometries),
indexes_species(multiphase[species]), K(K)
{
	SetKinetics(kinetics);
}

template<class KineticsType> void
Reaction::SetKinetics(const KineticsType& kinetics)
{
	this->kinetics = std::shared_ptr<Kinetics> (new KineticsType(kinetics));
}

ostream&
operator<<(ostream& out, const Reaction& reaction);

#endif /* REACTION_H_ */
